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Abstract 

In this introductory review, we give an overview of the computational chemistry methods commonly used in the field of metal-organic 
frameworks (MOFs), to describe or predict the structures themselves and characterize their various properties, either at the quantum 
chemical level or through classical molecular simulation. We discuss the methods for the prediction of crystal structures, geometrical 
properties and large-scale screening of hypothetical MOFs, as well as their thermal and mechanical properties. A separate section 
deals with the simulation of adsorption of fluids and fluid mixtures in MOFs. 
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Since its emergence in the 1950s, molecular simulation has 
seen an ever-growing use in research in the fields of physical, 
chemical, and materials sciences, where it offers an additional 
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year topic ref. 

2015 general modeling of MOFs (book) Ill 

2015 quantum chemical characterization of MOFs, including catalysis lO 

2014 high-throughput computational screening m 

2014 hrst-principles force helds for guest molecules @1 

2013 gas separation m 

2012 screening for adsorption and separation ||6l 

2012 methane, hydrogen, and acetylene storage (71 

2011 adsorption in flexible MOFs (H 

2011 energy, environmental and pharmaceutical applications (91 

2011 screening for separation applications da 

2009 hydrogen storage (TTl 

2009 adsorption (TTl 

2008 adsorption and transport (TSl 

2007 adsorption of small molecules (T4l 


Table 1: List of reviews published on computational characterization of metal-organic frameworks. 


dimension to the characterization and understanding of systems, 
complementary to experimental techniques and pen-and-paper 
theoretical models. The development of novel computational 
methodologies, together with the exponential increase in com¬ 
putational power available to researchers, have dramatically 
expanded the range of problems that can be addressed through 
modeling. This is true of resource-intensive calculations per¬ 
formed on high-performance computing (HPC) supercomputers, 
but it is also true of desktop workstations, and even now of 
laptops and mobile devices. csiiia Several of the simulation 
techniques of computational chemistry have now reached the 
status of being relatively “routine” calculations and are nowa¬ 
days considered an integral part of the researcher’s toolbox, just 
like experimental characterization techniques like X-ray diffrac¬ 
tion and NMR spectroscopy. Among those, we can cite Density 
Functional Theory calculations and Grand Canonical Monte 
Carlo simulations. It is possible to become a user of these tools 
with relatively little training, relying either on commercial or 
academic software with user-friendly interfaces. However, as 
with any technique, one should always take great care in check¬ 
ing the validity of the tools for the system at hand, as well as in 
interpreting the results obtained. 

Given the formidable research effort focused on metal- 
organic frameworks (MOFs) in the past decade, with more 
than 20,000 papers published (at a current rate of more than 
5 MOF papers per day (TTl), 15,000 structures on record at the 
Cambridge Crystallographic Data Centre, and over 170 review 
articles dedicated to this topic, the published literature on com¬ 
putational studies of MOFs is in itself abundant. Theoretical 
approaches are in many cases used, in combination with experi¬ 
mental characterization techniques, to study newly synthesized 
materials and understand their properties at the microscopic 
level. In this introductory review, we give an overview of the 
computational chemistry methods commonly used in the field 
of MOFs, to describe the structures themselves and character¬ 
ize (or predict) their various properties. It is by no means a 
systematic review of existing computational work on MOFs, of 
which there simply is too much to systematically enumerate here. 


Rather, we will try to give the reader an idea of what is possible 
in theoretical studies of MOFs, both at the quantum chemical 
level and through classical molecular simulations. For more 
details on specific areas of interest, we refer to existing revie ws 
of computational MOF studies, which are listed in Table 1 ^ We 


finish the review by pointing out some of the open questions and 
challenges in the field. 


2. Structural properties 

2.1. Crystal structure prediction 

In most cases, the structure of a newly synthesized metal- 
organic framework is determined experimentally, using single¬ 
crystal X-ray diffraction data when possible, or solving the 
structure from powder diffraction data if single crystals of suf¬ 
ficient size or quality cannot be obtained. In the latter case, 
it often happens that because of the molecular complexity of 
the material, a low symmetry, or a large unit cell, the structure 
solution is arduous or impossible to solve. Moreover, in other 
cases, there is a need for methodologies of true (or ab initio) 
computational prediction of MOF crystal structures, without any 
input of experimental data. These structures can then either be 
used to guide the synthesis of novel materials, or to identify 
materials already synthesized but whose crystals structure has 
not been solved experimentally, by comparison of X-ray powder 
diffraction patterns. 

In the past decades, computational crystal structure prediction 
has made giant strides in both the fields of molecular crystals (T^ 
and inorganic materials (dense and porous alike). im A large 
variety of methods have been developed and perfect to exploring 
the configurational space of these materials, including algo¬ 
rithms based on the simulated annealing approach, genetic (or 
evolutionary) algorithm methods. (201 ITTIl and molecular pack¬ 
ing approaches. (221 There have also been techniques developed 


Tn particular, we do not address in this review the very specific topic of the 
catalytic activity of MOFs: on this, we refer the reader to the very recent review 
by Odoh et al.(21 
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that are targeted specihcally at framework materials, including 
both extended inorganic solids (such as zeolites) and hybrid 
inorganic-organic frameworks. Indeed, framework materials 
present a specihc challenge when it comes to structure prediction, 
featuring both strong directional bonds and weaker dispersive 
interactions. We briefly summarize here the crystal structure 
prediction commonly used for metal-organic frameworks, and 
refer the reader to existing reviews (231 |24|[25l for a more com¬ 
prehensive treatment of the subject. 

Metal-organic frameworks, like inorganic framework mate¬ 
rials, can be regarded as composed of elementary or secondary 
building units (SBU) assembled together into three-dimensional 
networks. (261 IZTl This fact is exploited in the Automated As¬ 
sembly of Secondary Building Units (AASBU) method for 
structure prediction. First developed for the computational pre¬ 
diction of inorganic extended lattices, (28l the AASBU method 
uses predefined SBUs with tailored “sticky-atom” interactions 
potentials between them, allowing the coordination in corner-, 
edge- and face-sharing modes. From these elements, the SBUs 
auto-assemble into three-dimensional frameworks through series 
of simulated annealing and energy minimizations. This strat¬ 
egy is applicable to MOFs as well as inorganic frameworks, as 
was demonstrated by the “prediction” of experimentally-known 
structures, including HKUST-1 and MOF-5, as well as novel hy¬ 
pothetical frameworks. (291 The AASBU method, which directly 
simulates in silico the assembly of MOFs, is rather expensive 
computationally, as it requires series of minimizations and simu¬ 
lated annealings in order to converge structures. 

An alternative method is the direct construction of frameworks 
from SBUs by direct enumeration of possible topologies com¬ 
patible with the “building blocks” available. In this top-down 
approach, sometimes called the “decoration strategy”, the po¬ 
tentially infinite space of polymorphs of a given composition 
is explored by decorating a list of mathematical nets with the 
chosen SBUs, which act as vertices and edges of the net. llSTlI^ 
The nets themselves can be enumerated systematically from 
mathematical algorithms, given certain constraints on their fea¬ 
tures and complexity. (33]| They can also be taken from a list of 
experimentally known structures, such as the database of known 
zeolitic structures (for zeolitic nets) (341 or the broader Reticular 
Chemistry Structure Resource (RCSR) database. (3511361 This 
approach has been used on a large variety of system, including: 

• the identification of carboxylate-based MOFs with zeotype 
topologies, including the complex MIL-100 and MIL-101 
structures; (37l |38l 

• the enumeration of hypothetical (or “not-yet-synthesized”, 
if one is optimistic!) Zeolitic Imidazolate Frameworks 
(ZIF),(39l EOl or more recently of porous zinc cyanide 
polymorphs : (40l 

• hypothetical covalent organic frameworks fCOF). (4Tl 

The structures generated through this decoration process then 
need to be “relaxed” through energy minimization, relying on 
classical force field-based simulations or quantum chemistry 
calculations, in order to confirm their stability and evaluate their 
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Figure 2: Depiction of the linker replacement strategy for computational predic¬ 
tion of isoreticular MOF structures, exemplified on the MIL-88 family. Repro¬ 
duced from ref. m with permission from The Royal Society of Chemistry. 


relative energy (or formation enthalpy). The latter is important in 
determining the experimental feasibility of structures: even if a 
given MOF structure is a local minimum in energy, if its relative 
energy is too high compared to other polymorphs, it cannot be 
considered a good target for possible synthesis. This process is 
exemplified on Figure[2 in the case of hypothetical ZIFs.|[39j|30l 
Starting with a given composition (in this case, Zn^+ cations and 
unsubstituted imidazolate anions), zeolitic nets from the IZA 
database (34l are decorated with the metal centers (replacing the 
silicon atoms) and organic linkers (replacing the oxygen atoms). 
These “ideal” starting structures are energy-minimized through 
DFT calculations. The resulting structures can then be further 
studied, and their formation enthalpy (here, their energy relative 
to dense polymorph ZlF-zni) give insight into their experimental 
feasibility. In the case of ZIFs, a general correlation between 
enthalpy and density is found, but the relatively small dispersion 
of the enthalpy values between hypothetical and experimental 
energies hints that the currently hypothetical frameworks should 
be attainable by solvothermal synthesis. 

Another property of metal-organic frameworks which can be 
leveraged for the computational prediction of new structures is 
the principle of isoreticularity, i.e. the possibility of a family of 
MOFs to share the same topology and metal centers, with varia¬ 
tions in their organic linkers both in terms of length and function¬ 
alization. This is best exemplified by the famous IRMOF series 
of structures, (421 which all possess the same topology as the 
original MOF-5 (aka IRMOF-1) structure. (43l Based on a given 
structure for a parent MOF, it is possible to predict the struc¬ 
tures of possible isoreticular analogues. Gradually replacing 
the original linker with longer organic molecules with the same 
coordination modes, and performing quantum chemical energy 
minimizations on each individual structure. This computational 
strategy is of interest in proposing structures of expanded ver¬ 
sions of known structures. This was demonstrated in the case 
of the MIL-8 8 family of materials, a series of iron (III)-based 
and chromium(III)-based MOFs with linear dicarboxylic acid as 
linkers. There, the ligand replacement strategy was used for 
computational structure elucidation of the MIL-88B, MIL-88C, 
and MIL-88D compounds (see Figure [^. (441 Isoreticularity 
was also leveraged to design and screen new MOF-5 analogues 
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Figure 1: Experimental and hypothetical zeolitic imidazolate frameworks (ZIP) generated by decoration of zeolite topologies by Lewis et aL, Eol along with a plot of 
their relative energy (compared to dense ZIP of zni topology) versus density. Adapted from ref. do) with permission from The Royal Society of Chemistry. 


based on commercially available organic linkers, da and for the 
de novo synthesis, after computational prediction, of ultrahigh 
surface area MOF NU-100. 1461 

As a hnal example, a combination of this linker replacement 
and functionalization approach with the “decoration” strategy 
highlighted above was used by Gomez-Gualdron et al. flTl to 
generate 204 hypothetical MOF structures. In this reverse 
topological approach. 1411 the authors used as starting SBUs a 
zirconium-based metal center (Zr 604 ), present in widely-studied 
material UiO-66(Zr), and 48 organic linkers (12 ditopic and 36 
tetratopic building blocks). These building blocks were then 
combined in four possible topologies: fcu,ftw, sen and esq (see 
Figure [^. After computational identiheation of a top performer 
candidate for methane volumetric deliverable capacity, the pre¬ 
dicted MOF NU-800 was synthesized and its high predicted gas 
uptake properties were conhrmed by experimental high-pressure 
isotherm measurements over a large temperature and pressure 
range. 

2.2. Large-scale screening, enumeration of hypothetical struc¬ 
tures 

Although the different methods of crystal structure prediction 
described in the previous section each have different strengths 
and weaknesses, it is clear that they cannot necessarily scale to 
generate very large numbers of hypothetical structures. Yet, the 
constant search for new metal-organic frameworks has lead to 
an intense research effort to replace the current trial-and-error 
approach to MOF synthesis by enabling computationally-guided 
design of novel MOFs. There has thus been an effort, in the past 
few years, to enable larger-scale screening of potential MOF 
structures by generating larger number of hypothetical structures, 
and storing them into databases along with basic characterization 
information. This effort has been further spurred by being part 
of the broader Materials Genome Initiative, a 100 million dollar 
effort from the White House that aims to '‘discover, develop, and 
deploy new materials twice as fast'’ as the current methods. 1481 
l49l In this section, we detail some of the recent milestones in 


s 

Pigure 4: Scheme of the bottom-up approach for large-scale generation of hypo¬ 
thetical MOFs by Wilmer et al. do). Reproduced from ref. 0 with permission 
from The Royal Society of Chemistry. 


the generation of large-scale databases of hypothetical materials 
and their screening for specific applications. O 

Hypothetical databases of zeolite structures have been avail¬ 
able for many years, EB112 Enin El containing up to two- 
million unique feasible structures. In the MOF world, the first 
large-size database of hypothetical material was the pioneer¬ 
ing work of Wilmer et al. in 2Q12. ll50ll a decisive step towards 
large-scale high-throughput screening of metal-organic frame¬ 
works. The generation procedure followed by Wilmer (depicted 
on Figure is iterative, based on step-wise combination of 
predefined building blocks relying on geometric rules of at¬ 
tachment, pretty much like LEGO® bricks. All this information 
input into the generation procedure comes from the experimental 
literature: the building blocks are based on the reagents used in 
reported MOF syntheses, and the geometric rules for attaching 
the blocks together are determined by the crystallographic struc¬ 
tures of known compounds using this particular coordination. 
This approach is more powerful than the decoration or linker 
replacement strategies, because it is entirely general and does 
not restrict to known topologies. It is also less computation¬ 
ally demanding than the assembly approaches such as AASBU, 
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Figure 3: Scheme of the reverse topological approach for crystal structure prediction, as implemented by Gomez-Gualdron et al. l47l to generate 204 hypothetical 
(Zr604)-based MOF structures. Adapted with permission from ref. ED- Copyright 2014 American Chemical Society. 


because it does not necessitate the minimization of energy (or 
scoring function). It thus scales to large number of structures, 
as Wilmer demonstrated by generating 137,953 hypothetical 
MOFs from a library of 120 building blocks, with the constraint 
that each MOF could contain only one type of metal node and 
one or two types of organic linkers, along with a single type of 
functional group. 

A second example of the generation of hypothetical MOF 
structures comes from the family of zeolitic imidazolate frame¬ 
works (ZIF). Based on the extensive database of more than 
300,000 hypothetical zeolite structures by Deem, ||55l[56]| Lin 
et al. have used the decoration strategy to obtain a database 
of Zn(imidazolate )2 ZIF structures, replacing the zeolites’ sil¬ 
icon atoms by zinc cations and the bridging oxygen atoms by 
imidazolate linkers. Ell In addition to these two hypothetical 
MOF databases, two publicly available databases of MOFs have 
been constructed from experimental structures. The hrst was 
gathered and published by Goldsmith et al. llSSl Derived from the 
Cambridge Structural Database (CSD),|[59l[60l the database con¬ 
tains computation-ready MOF structures, i.e. disorder-free 
structures from which solvent has been removed. The authors 
originally used it for the screening and selection of optimal hy¬ 
drogen storage materials, as well as highlighting the theoretical 
limits of hydrogen storage in MOF materials. 

A second database was built recently by Chung et al., with 
the aim of producing ''a nearly comprehensive set of porous 
MOF structures that are derived directly from experimental 
data but are immediately suitable for molecular simulations or 
visualization”. len This database of computation-ready, experi¬ 
mental MOF structures (or CoRE MOF database) was produced 
from the CSD crystallographic structures through a multi-step 
procedure summarized in Figure (i) selection of potential 
MOFs, based on chemical bond analysis (over 60,000 candidates 
after this step); (ii) elimination of 1-D coordination polymers and 
2-D hydrogen bonded networks to retain only three-dimensional 


MOFs (over 20,000); (hi) removing partially occupied crystallo¬ 
graphic positions, elimination of structures with disorder in the 
framework, identification of charge-balancing ions; (iv) solvent 
removal. The final CoRE MOF database contains over 4,700 
hypothetical porous MOF structures. Its authors used it to inves¬ 
tigate the structural properties of the CoRE MOEs that govern 
methane storage capacity in MOFs. 

All the databases described above, whether of experimental 
and hypothetical structures, including zeolites, MOFs, porous 
polymer networks, and more, have been consolidated as part 
of The Materials Project and available online (upon registra¬ 
tion) for browsing as well as systematic data mining through 
documented Application Programming Interfaces (APIs). ll6^ 
Their use as a basis for large-scale high-throughput screen¬ 
ing of porous materials has taken off in the past three years. 
Most of the screening studies proposed so far focus on iden¬ 
tifying materials for adsorption-based applications, including 
separation. |[63l l64l capture. l(57l and storage ISOl 1^ [66l l67l of 
strategic gases: hydrogen, | 66 l| 62 l methane, EollMlISl carbon 
dioxide, lEil noble gases, (631 [64l etc. In these high-throughput 
screening approaches, the performance of materials are typically 
evaluated either by geometrical properties (pore size, pore vol¬ 
ume, accessible surface area; see Sectioi jZ3| or by evaluation of 
the host-guest interactions and adsorption properties (through 
Grand Canonical Monte Carlo simulations with classical force 
fields; see Section [4T] ). Large-scale high-throughput screen¬ 
ing of metal-organic frameworks targeting other properties, or 
relying on other descriptors, has not developed so far. 

2.3. Geometrical properties: surface area, pore volume and 
pore size distribution 

Once a MOF structure has been determined, either crystallo- 
graphically or using a combination of diffraction experiments 
and quantum chemistry calculations, there are several so-called 
geometrical properties that can be calculated based solely on 
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Figure 5: Schematic illustration of the CoRE (Computation-Ready Experimental) MOE database construction. (6T] Adapted with permission from ref. (eT). Copyright 
2014 American Chemical Society. 


the unit cell parameters and atomic positions. Such geometrical 
calculations are straightforward to perform and computationally 
inexpensive, ranging from seconds to minutes on a desktop com¬ 
puter. They rely on a description of the material where each atom 
of the framework is a hard sphere, centered on crystallographic 
coordinates. The radii used for these hard spheres are the van 
der Waals radii of the atoms. Probe molecules, used in order to 
calculate guest accessibility, are also considered spherical, with 
a diameter equal to the kinetic diameter of the molecule. 

The most common geometric characterization performed on 
microporous solids is that of their accessible surface area and 
accessible pore volume, representing the surface (resp. volume) 
accessible to guest molecules of a given size. Three possible 
surfaces can be used to measure this, as depicted in Figurej^ The 
van der Waals surface is that dehned by hard spheres centered on 
each atom, and does not depend on probe size. The accessible 
surface is that accessible to the center of mass of a probe of 
radius r^. Finally, the Connolly surface (also called solvent- 
excluded surface) delimits the space which is inaccessible to any 
part of the spherical probe; it is technically harder to compute. 
Diiren et al. have shown that the accessible surface area is the 
appropriate surface area to characterize crystalline solids for 
adsorption applications. 1691 In particular, the accessible surface 
area calculated with a probe size corresponding to nitrogen 
(rp = 3.681 A fTOl ) can be directly compared with BET surfaces 
from experimental nitrogen isotherms ITTl (though experimental 
values might be lower because of blocked pores, or higher in 
presence of defects; see Section [46| ). 

Sampling methods for the calculation of molecular surface 
are the most common numerical methods to evaluate the solvent- 
accessible surface of both molecules and periodic crystalline 
systems. This procedure is illustrated on Figure [^s lower panel: 
from a sample of points on each atom’s “solvation sphere” (of 


radius + rp), the proportion of points that are not buried inside 
neighboring atoms determines the atom’s contribution to the 
total accessible surface area. This is called the Shrake-Rupley 
algorithm. Cl A similar sampling method can be adopted for 
the accessible pore volume: the pore space is sampled, e.g. on 
a regular mesh, and the number of mesh points falling within 
the pore volume is counted. m Other sampling methods exist, 
however, for the determination of both accessible surface area 
and accessible pore volume, such as the use of ray casting |l74ll 
and analytical calculations. ITSl 1761 

Another tool at our disposal for geometric characterization 
of pore space is the pore size distribution (PSD), as illustrated 
on Figure in the case of metal-organic framework HKUST-1 
(also known as Cu3(btc)2). Experimentally, pore size distribu¬ 
tions can be obtained by numerical analysis of experimental 
low-temperature nitrogen or argon adsorption isotherms, mi eg 
given the choice of a reference pore geometry (slit-like, cylin¬ 
drical, spherical) and of an approximate chemical composition 
(though no kernels are available specifically for MOF materials). 
Thus, they provide a point of qualitative comparison between 
the geometric properties of the ideal crystal structure and those 
of the sample’s pore space as observed through the adsorption 
process. The standard method to calculate geometrical pore 
size distributions was developed for Vycor glasses by Gelb 
and Gubbins. llSOll and its computational efficiency was later im¬ 
proved by the same group. IMl This method is entirely generic, 
applicable to nanoporous materials containing both micropores 
(below 2 nm) and mesopores (between 2 nm and 50 nm) and has 
been used with success in several MOF materials. 1821 [8^ 

2.4. Advanced geometrical descriptors 

These methods described so far, however, do not account 
for the connectivity of the pore space determined by geomet- 
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Figure 7: Left: view on the structure of HKUST-1 (Cu 3 (btc) 2 ), featuring small pores (green) and two different types of large pores (green and red). Right: pore size 
distribution (PSD) calculated from the crystal structure of HKUST-1, with peaks corresponding to the three types of pores. Adapted from ref. (77], with permission 
from Elsevier. 




Figure 6: Upper panel: definition of the van der Waals surface (black line), 
Connolly surface (blue line), and accessible surface area (red dashed line). 
Lower panel: scheme of the calculation of accessible surface area contribution 
by each atom through sampling of its “solvation sphere”: green points are 
accessible, red points are buried inside neighboring atoms and thus inaccessible. 


ric means, nor for its accessibility to guest molecules along a 
diffusion path. For example, if a nanoporous structure presents 
a cavity linked to a single channel (side pocket), a guest may 
ht inside the cavity but not be able to enter in the hrst place 
if the channel is too narrow. It is therefore necessary to ana¬ 
lyze the pore volume (and associated surface area) by breaking 
it down into path wise connected components. Components 
of dimensionality equal to zero correspond to isolated cages 
or unreachable side pockets, which sorbate molecules cannot 
dynamically enter. Once identihed, these need to be excluded 
from GCMC simulations of adsorption, in order to avoid inser¬ 
tion of molecules where it is not physically realistic. ll84l [85l 


Components of higher dimensionality are one-dimensional pore 
channels, and 2D and 3D pore networks. 

The sampling methods described above for accessible surface 
and pore volume determination can be extended to further parti¬ 
tion the porous network into connected components, typically by 
mesh-based propagation methods lISTl [8^ [8^ [83]| derived from 
the classical algorithms developed in the study of percolation 
theory. (901 These grid-based methods are computationally ex¬ 
pensive, in particular for high-accuracy determinations which 
require very hne mesh spacing. An alternative to the use of grid- 
based sampling to characterize pore accessibility exists, based 
on Voronoi decomposition (depicted on Figure [^,|[9TJ|92l which 
for a given arrangement of atoms in a periodic domain provides 
a graph representation of the void space. The analysis of the 
Voronoi network and the accessibility of its nodes can then yield 
information into the components of the pore system, the dimen¬ 
sionality of the different channel systems, and the associated 
surface area and volume. IMIl 

In addition, the analysis of the Voronoi network can yield 
other useful geometric descriptors for porous materials. Three 
quantities, in particular, are of particular relevance to the descrip¬ 
tion of molecular adsorption and transport in porous materials 
(see Figure[g):ll93llF7ll^ 

• the diameter of the largest included sphere, Di, which 
reflects the size of the largest cavity within a porous mate¬ 
rial; 

• the diameter of the largest free sphere, Df, representing 
the largest spherical probe that can diffuse fully through a 
structure (i.e. the size of the narrowest constriction in the 
channel system); 

• the largest included sphere along the free sphere path, 

Af. 

The definitions of these three diameters are illustrated on Fig¬ 
ure They can be used in order to better understand adsorption, 
separation and diffusion of guests in nanoporous materials, 19^ 
to quantify the similarities (and differences) between pore spaces 
of materials in a given family, as well as for high-throughput 
screening of structure databases. |[86l[3]| 
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Figure 8: Advanced geometrical characterization of nanoporous materials. Top left: depiction of the Voronoi decomposition of space in 2D and 3D. Top right: 
definition of the diameter of the largest included sphere (A), diameter of the largest free sphere (Df), and largest included sphere along the free sphere path (Af)- 
Bottom: decomposition of the pore volume of the DDR zeolitic framework (for a probe of radius 3.2 A) between accessible and inaccessible volumes. Adapted from 
ref. f^ . with permission from Elsevier. 


From a practical point of view, the Voronoi-based analysis 
of nanoporous materials described above is implemented in the 
open-source Zeo-f-i- software package. 1951 [86l [961 It allows in 
one single run the calculation of all geometric features of a 
given system. It is widely used for large-scale studies of zeolite 
and MOF databases and the screening of hypothetical novel 
structures. (961 

Finally, while Voronoi-based descriptions of pore space ap¬ 
pear as the most generic tool for systematic description of pore 
space geometries, it is worth noting that other approaches have 
been developed. We will cite here in particular the alternative 
method of First et al., which aims at fragmenting the pore space 
and representing as a set of geometrical blocks such as cylinders 
and spheres in order to identify portals, channels, cages, and their 
connectivity. 1971 [98l Performing this analysis on a large number 
of available experimental crystal structures. First et al. built two 
online databases of nanoporous materials, ZEQMICS l99l and 
MQFQMICS. 1 1 001 aggregating quantitative information on the 
geometrical characteristics of zeolites and MOFs, respectively. 

2.5. Localization of extra-framework ions 

While most metal-organic frameworks possess a neutral 
framework, there is an important subclass of MOFs that feature 


ionic frameworks and charge-compensating extra-framework 
ions in their pores. These materials present some similarity 
to cationic zeolites, although in the case of MOFs both an¬ 
ionic frameworks (with extra-framework cations) and cationic 
frameworks (with extra-framework anions) are possible. In 
the past decade, a large number of ionic MOFs (sometimes 
called charged MOFs) syntheses have been reported in the 
literature. ill 04i Perhaps the best-known materials in this fam¬ 
ily are the anionic zeolite-like metal-organic frameworks, or 
ZMOFs. 11051 rho-7MO¥ and sod-7MO¥, synthesized by Ed- 
daoudi and coworkers, are porous anionic MOFs with ze¬ 
olitic topologies, whose overall neutrality is ensured by charge¬ 
balancing 1,3,4,6,7,8-hexahydro-2if-pyrimido[l ,2-a]pyrimidine 
cations. cna 

Both cationic and anionic metal-organic frameworks have 
been studied for potential applications. Due to the electric 
fields generated inside their nanopores by the presence of extra- 
frameworks ions, ionic MOFs show strong interactions with 
guest molecules, and specific interactions with polar guests in 
particular. This effect can in turn be leveraged for applications 
in gas adsorption and storage. ! 1071 1 1081 separation, 1 1091 molec¬ 
ular recognition, nioi drug delivery. !! 11! ion exchange !! 12! 
and catalysis. !! 1311114111 15! Since the extra-framework ions in 
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Figure 9: Analysis of the pore systems of zeolite MFI and metal-organic 
frameworks NOTT-401. fToTl ZIF-7. IT02l and Mg-r/io-ZMOF lT03l by the 
ZEOMICS(92l and MOFOMICSl^ tools. Reproduced from ref. |92l with 
permission from The Royal Society of Chemistry, and the MOFOMICS 
website. ITool 
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Figure 10: Top: localization of the two cationic sites (I and II) for Na"^ in the 
anionic rho-ZMOF framework. Middle: Coadsorption isotherms and selectivity 
for an equimolar CO 2 /CH 4 mixture, calculated by Grand Canonical Monte Carlo. 
Bottom: Locations of CO 2 molecules adsorbed from the equimolar CO 2 /CH 4 
mixture, at various pressures. Reproduced with permission from ref. IT22I . 
Copyright 2009 American Chemical Society. 


charged MOFs can be exchanged, these materials offer a large 
versatility and tunability. However, the rational design of novel 
materials and their post-synthetic optimization requires a good 
understanding of their behavior at the microscopic level, and in 
particular of the localization of their extra-framework cations 
and the nature and strength of the ion-guest interactions. 

The localization of extra-framework ions in the charged MOF 
structures is not always possible, in particular because the ions 
may be too delocalized or disordered. There, molecular simula¬ 
tion can play an important role, identifying possible binding sites 
for ions (i.e. local energy minima) and their physical character¬ 
istics: binding energy, mobility (i.e. spatial spread of the site). 
This can be achieved either at the level of quantum chemical cal¬ 
culations, or through molecular simulations based on classical 
force fields. The latter, in addition to individual binding sites, can 
also help determine the distribution of cations between the vari¬ 
ous possible sites as a function of the experimental conditions: 
number and nature of cations, temperature, presence of solvent, 
etc. This is typically achieved through Monte Carlo simulations 
including large-scale displacement moves (or “jumps”), in order 
to overcome the formidable energy barriers typically involved 
with movement of a cation from site to site. These simulations 
can also borrow methodologies from the very extensive scien¬ 
tific literature available on the topic of extra-framework cation 
localization in zeolites. Hl 16111171 [11811 including the use of sim¬ 
ulated annealing 1 1191 or parallel tempering methods !! 1201112 111 


to reach equilibrium in reasonable time. 

Once the localization of extra-framework ions has been de¬ 
termined, it can then be used as input or basis for the study 
of adsorption properties of the ionic MOF. 1 1231 II 241 [T25l 1 1091 
These are methodologically similar to simulations of adsorp¬ 
tion in neutral MOFs, except for the very large strength of the 
Coulombic ion-guest interactions. While simulations of adsorp¬ 
tion are treated in detail in Section we give in Figure a 
rather typical example of results from molecular simulation in 
ionic MOFs, here in the case of the anionic rho-ZMOF with 
extra-framework Na"^ cations. In this study, Babarao et al. II 12211 
identified two types of binding sites for Na"^ ions in the anionic 
framework, with site I in an eight-membered ring and site II in 
the a-cage. The authors showed that carbon dioxide is adsorbed 
predominantly over other gases, including methane, because of 
its strong electrostatic interactions with the charged framework 
and the presence of Na"^ ions acting as additional adsorption 
sites. 

3. Physical properties 

3.1. Mechanical properties 

Compared to structural characterization and adsorption prop¬ 
erties, studies of the mechanical properties of metal-organic 
































Published as: Coord. Chem. Rev. 2015, DOI: 10.1016/j.ccr.2015.08.001 


14 

12 

10 

8 

6 

4 

2 

0 



Minimal Shear 
Modulus (GPa) 











□ □ □ 



MOF-5 



chanical stability. The UiO -66 family of materials is one of the 
exceptions, with shear modulus in the 12-14 GPa range, due to 
its strong Zr-0 linkages and high degree of coordination. II 1 2611 
Characterization of elastic constants of MOFs can also be per¬ 
formed on crystal structures without cubic symmetry, with the 
same techniques. Although the manual generation of strained 
structures by hand makes it more tedious in lower-symmetry 
crystal classes, some quantum chemistry software now handle 
it directly (including CRYSTAL14 II140[ [14111 ). and scripts are 
available for others. These allow to calculate the full four- 
rank tensor C of second-order elastic constants, which re¬ 
lates strain £ to stress G in the tensorial Hooke’s law: 

^ij = ( 1 ) 

kl 

The number of independent nonzero elements of this stiffness 
tensor, depends on the crystal class. It determines the entire 
elastic behavior of the material, and by tensorial analysis can be 
used to calculate physical properties of interest (see Figure [T^, 
including: 


Figure 11: Top: comparison of the minimal shear modulus of several metal- 
organic frameworks, obtained from DFT calculations. Bottom: representation 
of the soft shear modes of MOF-5 and HKUST-1. Reproduced with permission 
from ref. Copyright 2013 American Chemical Society. 


frameworks appeared relatively late in the literature, both ex¬ 
perimental and computational. The hrst theoretical calculations 
of mechanical properties were predictions of the bulk modulus 
of MOF-5 (and several analogues) by DFT calculations: 1 1271 
Fuentes-Cabrera et al. calculated the energy vs. volume curves 
for each of these cubic materials, fully relaxing the atomic posi¬ 
tions for each value of unit cell parameter a\ then the curves were 
htted by the Birch-Murnaghan equation of state. II 12811 Later 
work performed calculations not only of bulk modulus, but also 
the individual elastic constants (Cn, C 12 , and C 44 ) of the cubic 
MOF-5. 1 1291 1 1301 [1311113211 along with its average Young’s and 
shear moduli. These calculations were again performed at the 
DFT level (either with LDA or GGA exchange-correlation func¬ 
tions). They employed either the htting of quadratic “energy vs. 
strain” curves, or linear “stress vs. strain” curves, by applying 
small strains to the relaxed structure. 

Since these seminal studies, other cubic materials have been 
studied by the same methods, including IRMOFs lll33L ZIF- 
8 II134L Ui0-66. II126II etc. In addition to these zero Kelvin quan¬ 
tum chemical calculations, other works have reported bulk mod¬ 
uli at finite temperature, through force field-based molecular dy¬ 
namics simulations, e.g. for MOF-5 1 1351 and HKUST-1. 1 1361 
The main conclusion, from both experimental and computational 
work on elastic properties of MOFs. 11137111381113911134II is that 
MOFs are far softer than inorganic nanoporous materials such as 
zeolites, i.e. they present much lower elastic moduli, though their 
detailed mechanical properties depend on both chenfical compo¬ 
sition and the framework’s geometric properties (see Figure pT]). 
In particular, many highly porous MOFs feature very low shear 
moduli (of the order of 1 GPa), implying relatively small me- 


• the directional Young’s modulus £’(u), also known as the 
tensile modulus, quantifies the deformation of the material 
in direction u, when it is compressed in that same direction; 

• the linear compressibility j3(u), which characterizes the 
compression along axis u when the crystal undergoes an 
isotropic compression; 

• the shear modulus (or modulus of rigidity) G(u, v), which 
quantifies the material’s response to shearing strains along 
u, in the plane normal to v; 

• the Poisson’s ratio v(u,v) which characterizes the trans¬ 
verse strain (in the v direction) under uniaxial stress (in the 
u direction). 

Programs are available for the calculation of these properties 
from the stiffness tensor, such as the ElAM II 14211 code or the 
online ELATE web app. ll43l 

The detailed analysis of the elastic properties of metal-organic 
frameworks has, in the last few years, been applied to different 
properties. In addition to the low elastic moduli of MOEs in 
general, Ortiz et al. 1 1451 [1431 demonstrated that the existence of 
high anisotropy in elastic properties, coupled with directions of 
very small Young’s and shear moduli (sub-GPa), is a signature 
of the flexibility of soft porous crystals A\A1\ determining their 
ability to undergo large structural deformations under stimula¬ 
tion (see Figure[T3]). 11481 This has allowed the prediction of flex¬ 
ibility for new structures, such as NOTT-300 and CAU-13 I149I 
(the latter has since been confirmed experimentally II 1 501 ). 

Among the mechanical properties of MOFs that have attracted 
some interest. Negative Linear Compressibility (NLC) I151I 
has drawn significant attention. This property of materials ex¬ 
panding in one or two directions under hydrostatic compres¬ 
sion is considered rather exotic in inorganic solids, lfT52l but 
has been observed in a relatively large number of framework 
materials ll53lll54lfT55lfT5^ and MOFs. fT57lfT58lfT5^ On the 
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Figure 12: Determination of the four-rank tensor C of second-order elastic constants, by quantum chemistry calculations, and its analysis to obtain physical properties 
such as Young’s modulus (E), shear modulus (G), linear compressibility (j8), and Poisson’s ratio (v). 


computational side, it can be studied in two ways. In the lin¬ 
ear elasticity regime, it can be studied through the computation 
of the elastic stiffness tensor as described above, as was done 
for the prediction and subsequent experimental conhrmation of 
negative linear compressibility in the MIL-53 family. II 15 81 It 
can also be performed by in silico compression experiments, 
studying the influence of finite increments of pressure on a 
MOF structure, either through enthalpy minimization calcula¬ 
tions under pressure, ihsti or through constant-pressure constant- 
temperature (A^, (T, r) molecular dynamics studies. 1 1491 

Because computational compression experiments allow the 
determination of structure (unit cell parameters and atomic po¬ 
sitions) as a function of applied mechanical pressure, it yields 
information outside of the elastic regime and can provide in¬ 
sight into the occurrence of pressure-induced structural tran¬ 
sitions. Two rather typical examples of these in silico compres¬ 
sions are shown in Figure [T^ using DFT-based enthalpy mini¬ 
mization of the crystal structures at increasing (or decreasing) 
values of pressure. The first is that of CAU-13. 11611 demon¬ 
strating linear elastic behavior at positive pressure and the exis¬ 
tence of a narrow-pore-to-large-pore structural transition under 
tension, i.e. for hydrostatic pressures around P —500 MPa. 
Although experimentally applying hydrostatic tension is not an 
option, it does correspond to the outward stress induced by ad¬ 
sorption of bulky molecules (such as xylenes) in small nanopores. 
The second example is that of ZAG-4, a Zinc Alkyl Gate mate¬ 
rial showing nonmonotonic behavior under compression, as seen 
from high-pressure single-crystal X-ray crystallography. csa 
Quantum chemical calculations of the compression process 
showed that the nonlinear behavior is associated with a struc¬ 
tural transition, namely a reversible pressure-induced proton 
transfer between an included water molecule and the linker’s 
phosphonate group. csni 

Finally, it should be noted that in the (relatively rare) case of 
materials for which a good flexible force field is available, such 
studies can also be performed using force field-based molecular 
dynamics simulations. This approach has been well demon¬ 


strated in the case of the crystal-to-crystal “breathing” transition 
in materials of the MIL-53 family (including MIL-47), as is dis¬ 
cussed fully in Section [43| The same approach was used to study 
the pressure-induced amorphization in ZIF-8 and ZIF-4. ifTMl 
There, fully-anisotropic constant-pressure molecular dynamics 
simulations were used to calculate the evolution of elastic con¬ 
stants of the materials at room temperature, as a function of 
pressure]^ This was used to show that the pressure-induced 
amorphization of both ZIF-8 and ZIF-4 is due to a shear mode 
softening under pressure. II 16311 leading to mechanical instabil¬ 
ity at sub-GPa pressures when the Bom stability conditions are 
no longer satisfied. II 16511 This approach was later extended to 
a larger database of ZIF structures, showing the very limited 
mechanical stability of a majority of both hypothetical and ex¬ 
perimental structures upon solvent or guest evacuation. II 16611 

3.2. Thermal properties 

Because of its scalar nature and limited practical range, the 
diversity in MOF response to temperature is somewhat nar¬ 
rower than their responses to pressure, described above. Yet 
there has been a relatively large number of computational 
studies on the thermal properties of MOFs, and in particular 
their thermal expansion. There is, among the whole class of 
metal-organic frameworks, a prevalence of negative thermal 
expansion (NTE) and its occurrence in relatively large tem¬ 
perature ranges often including room conditions. Materials 
exhibiting NTE contract when heated, a rare property among 
dense inorganic materials and usually limited to certain types of 
stmctures. IEtI Negative thermal expansion, however, is quite 
common among molecular frameworks, and has been observed 


^The elastic constants Qj can be obtained from the fluctuations of the unit 
cell vectors in constant-stress {N, cr, T) molecular dynamics simulations with 
anisotropic variations of the unit cell, using the strain-fluctuation formula: fTeil 
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Figure 13: Anisotropy of the elastic moduli as a signature of soft porous crystals: 
3D representation of the directional Young’s modulus of “breathing” metal- 
organic framework MIL-47, compared to MOF-5. On the bottom panel are 
indicated the stiffest and softest directions of deformation (minimal and maximal 
Young’s modulus). Adapted with permission from ref. 11441 . Copyright 2013 
American Chemical Society. 


experimentally in many metal-organic frameworks, including 
HKUST-l.lfT^ MOF-5, (163 other members of the IRMOF 
family. II 17 OB and many ZIFs. lll66B 

One possible way to study the thermal expansion of MOFs 
is to perform constant-pressure constant-temperature (N, cr, T) 
molecular dynamics simulations for various values of tempera¬ 
ture. The evolution of the volume and unit cell parameters then 
allow a direct determination of both the volumetric and linear 
thermal expansion coefficients, ay and a£ respectively: 



Figure 14: Examples of in silico compression experiments, by DFT-based 
enthalpy minimization under pressure. Top: narrow-to-large-pore structural 
transition in CAU-13 under tension. Bottom: pressure-induced proton jump in 
ZAG-4. Adapted with permission from refs. EH (from The Royal Society of 
Chemistry) and ITeol (Copyright 2014 American Chemical Society.) 


their thermal expansion is to calculate the thermal expansion 
coefficients within the quasi-harmonic approximation. This 
requires the determination of the phonon modes and frequencies 
of the structure, and their dependence on unit cell parameters: 
this is typically done by ab initio lattice dynamics calculations of 
the phonon frequencies at various points throughout the Brillouin 
zone. For each phonon mode i and k point, a mode-specihc 
Griineisen parameter, can be calculated: 


1 fdV\ ^ f di\ 


y 




( 2 ) 


YiM 


( ^log{0,,k \ 
V dlogV ) 


(3) 


where ^ is any unit cell parameter. This was used to conhrm 
the existence of NTE in several materials of the IRMOF family 
(using force held-based MD simulations), as well as to shed 
light onto its microscopic origin. lll35lll70B The propensity of 
IRMOFs for NTE was attributed to the presence of many soft 
(low frequency) transverse vibrational modes in the frameworks, 
and in particular the vibration modes of the linkers transverse 
to the metal-metal axes, resulting in a shorter average effective 
length of the linkers, even though all the bond lengths in the 
structure increase with temperature. In a later study, Peterson et 
al. combined neutron scattering with ab initio molecular dynam¬ 
ics of the metal clusters (copper paddle wheels) of HKUST-1 to 
elucidate the origin of its NTE, involving both concerted trans¬ 
verse vibrations as well as local molecular vibrations such as 
paddle wheel twisting deformation. Clll 

Another approach to study the lattice dynamics of MOFs and 


The phonon modes with both negative Griineisen parameter and 
low frequency (and thus large amplitude) drive the negative ther¬ 
mal expansion behavior. This approach was used, for example, 
to study the NTE in MOF-5, first at the F-point onlv. ll69B then 
later across the full Brillouin zone. iml This allowed to identify 
both optic and accoustic modes contributing to the macroscopic 
NTE of MOF-5 (as schematized in Figure p~5]). Unlike the direct 
molecular dynamics approach, the quasiharmonic approxima¬ 
tion cannot account for high anharmonicity, but it does give 
more insight into the roots of the negative thermal expansion at 
the microscopic level. 

Finally, another thermal property of interest in MOFs is their 
heat capacity, an important thermodynamic parameter to char¬ 
acterize the material itself and also optimized the adsorption 
process for practical applications. Porous solids with larger heat 
capacity can better adsorb the heat generated during adsorp- 
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Figure 15: Top: MOF-5 unit cell and main NTE-contributing vibration modes of the linkers. Bottom: phonon dispersion curves and density of states (bottom left), and 
associated values of Griineisen parameters (bottom right; redder = more negative). Adapted from ref. IT721 with permission from The Royal Society of Chemistry. 


tion, thus limiting the undesirable thermal effects (increased 
temperature leading to lower adsorption). However, very few 
papers directly address this issue, whether experimentally H 173 II 
or computationally. lEl 

3.3. Optical and electronic properties 

Optical properties of metal-organic frameworks, including 
UV-visible absorption and luminescence, have attracted a great 
deal of attention due to their potential for applications in chem¬ 
ical, biological, and radiation detection, medical imaging, and 
electro-optical devices. HI761 1 17711 In addition, there is also great 
interest in the characterization and control of band gaps, for 
applications in electronic devices, solar energy harvesting and 
photocatalvsis. I1781I179II180]| These optical and electronic prop¬ 
erties of MOFs, require modeling at the quantum chemical level, 
mostly in the Density Functional Theory (DFT) approach, in 
order to fully characterize the electronic structure of the material. 

A lot of the early work on the modeling of electronic 
properties of MOFs focused on the characterization of their 
band gaps J129l and the possibilities for tuning these by 
metal exchange (127] I181II or ligand functionalization H 1 8211 or 
substitution. 1 1831 Probably one of the most comprehensive ex¬ 
amples of this kind of studies is that of Hendon et ah 1 1841 
who explored through a combination of synthetic and computa¬ 
tional work the influence of linker functionalization on the band 
gap of MIL-125, a photochromic MOF based on Ti02 and 1,4- 
benzenedicarboxvlate. H 1 8511 By studying materials built from 
linkers with various functional groups, the authors demonstrated 
that the diaminated linker bdc-(NH 2)2 was the most efficient in 
lowering the band gap, from 3.6 to 1.3 eV. They also demon¬ 


strated that the introduction of just one aminated linker per unit 
cell was sufficient in lowering the band gap. 

In addition to the value of band gap, more recent work has 
shifted focus to a broader range of electronic properties, in¬ 
cluding detailed analysis of electron density and electrostatic 
potential. One such example is the recently proposed method 
by Butler et al. to report vertical ionization energy of MOFs, 
with respect to a vacuum level set at the value of the electrostatic 
potential (for MOFs whose pores are large enough). Ea Based 
on band gap and ionization potential values, the authors explain 
certain electrochemical, optical, and electrical properties of the 
materials studied. This method was also used for the charac¬ 
terization of the piezochromism of MOFs, i.e. their ability to 
change electronic and optical properties as a function of applied 
pressure. 1118611187II 

Finally, another area of interest is that of the luminescence of 
MOFs. The experimental literature on the topic (see for example 
the reviews in Refs. 117611177l and ll88b largely outweighs the 
computational studies, and this is an area still very much in 
development. Several studies have attempted to determined the 
electronic nature of the absorption and emission transitions in 
MOFs, and their dependence on linker functionalization, nature 
of the metal center, and linker-linker interactions and stacking. 
The accurate determination of the luminescence properties re¬ 
quire the use of computational intensive time-dependent density 
functional theory (TDDFT) H 189 11 1 90l with ad hoc exchange- 
correlation functionals and large basis sets in order to describe 
the excited states and optical (absorption and emission) spectra 
of MOFs. This has allowed to determine the nature of the emis¬ 
sion transition in a few luminescent MOFs J191lll92l including 
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Figure 16: Left: calculation of MIL-125’s electrostatic reference potential for MOFs, depicting the reference point chosen at the pore center. Right: vertical ionization 
energy of six porous MOFs with respect to a common vacuum level. Reproduced with permission from ref. ED Copyright 2014 American Chemical Society. 


archetypical MQF-5. il 1931 4. Adsorption 

As nanoporous materials, adsorption of molecular fluids in¬ 
side MOFs is one of their earliest and most studied properties. 
Due to their well-defined crystalline structure, high pore vol¬ 
ume, large surface area and tailorable pore sizes, MOFs show 
great promise in becoming the next generation of nanoporous 
adsorbents. They are thus natural candidates to supplement or 
replace zeolites in adsorption-based industrial applications, with 
a specific focus on energy-related and biomedical applications. 
This includes drug delivery, gas adsorption and capture, gas stor¬ 
age and delivery, separation in gas and liquid phase, purification, 
and sensing. 1 1941 1 195 1 11961 In particular, the separation, capture, 
and sequestration of carbon dioxide from industrial and auto¬ 
motive emissions has recently attracted intense research interest 
in the effort to curb emissions of this greenhouse gas linked to 
human-induced global warming. 1 197111981 

Adsorption properties of new porous materials are routinely re¬ 
ported along with the synthesis, structure and other physical and 
chemical characterization. In addition to the standard nitrogen or 
argon adsorption at cryogenic temperatures, part of the BET mea¬ 
surement of surface area J1991l200ll adsorption and desorption 
isotherms of gases of strategic interest such as CO 2 , CH 4 , CO, 
N 2 , and O 2 , is often performed to assess novel materials. Some¬ 
what more specialized topic, which are nonetheless of practical 
importance, such as adsorption of hydrogen. ifTTl 1^111202112031 
polar molecules (such as water 1120411205 1 and alcohols 1120611 ). 
and larger molecules (including hydrocarbons 12071 ) have all 
been extensively studied in some of the materials. 

Given the importance of this field, there has been significant 
computational work addressing both the fundamental under¬ 
standing of adsorption in metal-organic frameworks and the 
practical applications. Molecular simulation of adsorption al¬ 
lows to shed light into the microscopic root of the behaviors 
observed experimentally. It can also be a powerful tool for the 
quantitative prediction of adsorption and coadsorption in MOFs, 
as well as the optimization of materials for adsorption-based 
processes. Reviews on this topic are refs. and Q. Here, 
we summarize the state of the art in molecular simulation of ad¬ 
sorption in metal-organic frameworks, starting with the classical 
Grand Canonical Monte Carlo methods for adsorption and coad¬ 
sorption of fluids, before listing advanced techniques developed 
in the last few years for the study of specific aspects of MOF 
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Figure 17: Schematic diagram of the Monte Carlo simulation of adsorption in a porous host, in the Grand Canonical thermodynamic ensemble (in orange) and in the 
osmotic ensemble (allowing deformation of the material, in green). 


adsorption. 

4.1. Grand Canonical Monte Carlo 

The standard molecular simulation technique to study the ther¬ 
modynamics of adsorption in rigid nanoporous materials is to 
use Monte Carlo simulations in the Grand Canonical ensemble, 
or Grand Canonical Monte Carlo (GCMC).|[208l|209l[2l0l It 
has been extensively used and validated in a large variety of 
nanoporous materials, including zeolites, nanoporous carbons, 
mesoporous materials, etc. This approach, schematized in Fig¬ 
ure [T^ models the adsorption of molecular fluids or mixtures 
inside a rigid porous matrix (system of fixed volume V), for 
given values of the temperature (T) and chemical potential of 
the fluid (/i). 

In Monte Carlo simulations, series of configurations of the 
system under study are generated in a stochastic process, by 
random moves accepted according to each configuration’s Boltz¬ 
mann probability. These molecular Monte Carlo moves typically 
include translation and rotation of molecules, as well as in¬ 
tramolecular displacements, i.e. changes in a given molecule’s 
conformation. Grand Canonical Monte Carlo simulations go 
beyond this by simulating the exchange of molecules between 
the interior of a material’s pore space and an external reservoir 
of bulk fluid, with additional moves consisting of insertion of 
molecules into the pore volume, as well as deletion (moving 
them back to the reservoir). This allows the direct simulation of 
an open thermodynamic ensemble, with thermodynamic equi¬ 
librium between two phases (bulk and adsorbate) without an 
explicit interface. 

The conditions of a Grand Canonical Monte Carlo simulation 
are close to the thermodynamic conditions during experimental 
adsorption measurements. Each point of an adsorption isotherm 
in GCMC is the result of a simulation at fixed (/i,y, T), and 
the full isotherm A^ads(M) is obtained by running simulations at 


different values of jl. This is close to experimental isotherms, 
which are typically measured as A^excess(^). where P is the pres¬ 
sure of the external fluid. Two differences between simulated 
and experimental adsorption isotherms need to be taken into 
account for quantitative comparison. The first is to relate the 
chemical potential jl (used in GCMC) to the pressure P (mea¬ 
sured experimentally) of the bulk fluid. This can be achieved 
through an equation of state, as {dil/dP)T =Vm{P^T), the mo¬ 
lar volume of the fluid. This equation of state for the fluid needs 
to be obtained from prior Monte Carlo simulations, or alterna¬ 
tively obtained from experimental data. In the specific case of 
low-pressure adsorption, the ideal gas law might be used, then 
equating pressure and fugacity: jii = jii^ RT ln{P/P^). The sec¬ 
ond important difference between isotherms calculated through 
GCMC and measured experimentally is that between absolute 
and excess adsorbed properties, respectively. This can be sig¬ 
nificant in high-pressure experiment, and needs to be accounted 
for before any comparison between theoretical and experimental 

data.Eni 

Finally, this short introduction to Grand Canonical Monte 
Carlo can be concluded with the simple question: what can 
one compute with GCMC simulations? The output of GCMC 
includes both equilibrium thermodynamic quantities, as well as 
a representative set of configurations of the system in the given 
conditions. The macroscopic quantities include the absolute 
adsorption uptake (A/ads), i-e. the average number of adsorbed 
molecules in the porous system, which is used in plotting ad¬ 
sorption isotherms. It also includes energetic quantities, such 
as the isosteric heat of adsorption (^st)- These quantities can 
be directly compared to experimental data for validation, or 
used predictively as input for thermodynamic models of indus¬ 
trial adsorption-based processes, such as fixed-bed adsorption 
columns. In addition, Monte Carlo simulations also result in the 
generation of a sample of representative configurations of the 
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system, from which structural information can be obtained. In 
adsorption, in particular, one can thus plot density distributions 
for adsorbed molecules, yielding microscopic insight into the 
adsorption mechanism. 

4.2. Classical interaction potentials 

Monte Carlo simulations, as described above, rely on the 
evaluation of the energy of each generated conhguration of the 
system in order to evaluate its Boltzmann probability and accept 
or reject it. In short, the accuracy of the averages computed 
from the GCMC are directly dictated by the accuracy of the de¬ 
scription used for the interactions of the molecules in the system. 
Because such energy calculations of single conhgurations will 
be made on the order of millions of times per GCMC simulation, 
quantum chemical calculation of the energy of each conhgura¬ 
tion is out of the picture. Thus, the interactions between the huid 
molecules and MOF need to be described using classical mod¬ 
els: interaction potentials, also called/orc^ fields, are analytical 
functions of interatomic distances. In this approximation, we 
typically break down the host-guest and guest-guest interactions 
into terms with analytical expressions and simple physical mean¬ 
ing. At the intermolecular level, these terms include Coulombic 
interactions, long-range dispersion, short-range interatomic re¬ 
pulsion, polarizability, etc. A wide variety of functional forms 
are available to describe these intermolecular interactions, in¬ 
cluding the common Lennard-Jones and Buckingham potentials, 
as well as the Morse potential and the Feynman-Hibbs quantum 
effective potential, necessary for including quantum effects in ad¬ 
sorption of light gases such as hydrogen. Each “force center”, on 
which these intermolecular potentials act, need not necessarily 
be an atom. Approaches such as United Atom tUA). 1121211 and 
its refinement Anistropic United Atom (AUA) i213l have shown 
that it is possible of grouping a functional group (CH, CH 2 , CH 3 , 
...) into a single force center. In addition to the intermolecular 
terms, intramolecular terms may need to be considered for all 
but the smallest, rigid molecules. Intramolecular terms typically 
include stretching, bending and torsion potentials. 

The choice of force fields is one between accuracy and trans¬ 
ferability. On the one hand, designing ad hoc force fields for 
every material and guest molecule studied is a very time consum¬ 
ing trial-and-error process, which some consider close enough 
to a black magic, but it can provide very accurate description of 
the interactions in the system. On the other hand, the so-called 
transferable (or universal) force fields, which describe the same 
types of atoms in many different materials with identical param¬ 
eters, are much simpler to use and adapt to new MOF structures 
but their simplicity comes at the cost of a more approximate 
description of the potential energy surface. 

When it comes to transferable force fields, the most used are 
UFF (Universal Force Field) Ell and DREIDINGEH. Both 
force fields are generic enough to provide potential parameters 
for elements throughout the periodic table, covering both the 
organic linkers and metal atoms of MOEs. Em Some groups 
have nevertheless worked on extending these universal force 
fields to handle more transition metals, or to improve the de¬ 
scription of certain types of atoms (like the oxygen atoms in 


metal oxide clusters). Enl In addition, when modeling the ad¬ 
sorption of polar guest molecules, these force fields need to 
be supplemented with atomic partial charges to describe the 
Coulombic interactions. Classically, several methods are avail¬ 
able for the calculation of partial charges in molecular systems 
based on high-level quantum chemical calculations, grouped 
broadly in two families: population analysis or charge parti¬ 
tioning methods (Mulliken. 1121811 Hirshfeld. 112191 Bader ll22QII ). 
and electrostatic potential fitting (such as CHelpG II221ll and 
RESP II222ll k Electrostatic potential fitting methods are inher¬ 
ently more suited for the purpose of obtaining partial atomic 
charges for interatomic Coulombic potentials, but were long 
limited to nonperiodic systems. The past five years have seen 
a lot of developments of the state of the art in this area, with 
a new generation of electrostatic potential-based methods for 
periodic solids such as REPEAT (Repeating Electrostatic Po¬ 
tential Extracted ATomic) 12231 and DDEC (Density-Derived 
Electrostatic and Chemical charges). II224II as well as their later 
refinements. 1122511226II 

Other approaches to the problem of atomic partial charges 
determination have been proposed. If one wants to perform rapid 
calculations at the detriment of quality of the results, classical 
approximations can be used in the form of charge equilibra¬ 
tion methods. tSOl I227II These are particularly useful for high- 
throughput characterization of large numbers of materials. At 
the other hand of the spectrum, it is also worth noting that the 
determination of atomic charges can also be bypassed entirely, 
using the full electrostatic potential map in the material’s unit 
cell, as determined from DET calculations. E2H1 This latter ap¬ 
proach can become, however, computational quite expensive for 
materials with large unit cells. 

While transferable force fields can lead to reasonable agree¬ 
ment in the case of adsorption of small nonpolar molecules in 
some MOEs, it should be noted again that they are a rather crude 
“first order” approximation of the interatomic interactions. Thus, 
there is a real need for a systematic methodology for developing 
reliable force fields for novel or hypothetical MOE structures. 
A large number of such methods have been proposed in the 
past few years, using computationally-demanding high-level 
quantum chemistry calculations as a basis for determining force 
field parameters. The case of partial charges has been discussed 
already above, but other terms can also be derived from quantum 
chemical calculations: short-range repulsion, dispersion, polar¬ 
ization interactions, etc. These terms are typically optimized 
by choosing a functional form beforehand, and optimizing all 
available force field parameters to match certain properties of 
the high-level calculations. Some methodologies focus on prop¬ 
erties of the equilibrium structure, such as lattice parameters, 
atomic positions, forces (i.e. first derivatives of the energy with 
respect to atomic coordinates), vibration frequencies (i.e. second 
derivatives), bulk modulus (second derivative of the energy with 
respect to strain), etc. Other methodologies rely on fitting energy, 
and sometimes atomic forces, on a large sample of representative 
configurations of the system, themselves extracted from prior 
molecular simulations. This sampling can be performed with 
MD or GCMC simulation using universal force fields, or ab 
initio molecular dynamics. 
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In conclusion, a lot of methodologies have been proposed for 
the development of MOF force helds from hrst principles, many 
with great success on limited number of materials. Some of the 
recent work in this area includes: MOF-FF J229I BTW-FF J174I 
QuickFK 112301123 1 II among others. 112321123 3 II However, none of 
them so far have demonstrated a decisive edge in universality 
or large-scale predictive power. Development of hrst-principles 
force helds remains for the most part an open challenge in the 
molecular simulation of MOF adsorption. More studies of sys¬ 
tematic evaluation of force held performance are needed. II234II 
and other options for parameter determination (such as genetic 
algorithms 1235 l[T36ll ) might be necessary to reach a truly uni¬ 
versal methodology. 

4.3. Coadsorption and separation 

The Grand Canonical Monte Carlo scheme, as described 
above, can model not only the adsorption of pure components, 
but also mixtures of molecular huids, by specifying the chemical 
potential of each component in the mixture: (/ii, /i 2 , • • •, F, T). 
Like simulation of pure phases, Monte Carlo modeling of mix¬ 
ture coadsorption can be performed to provide microscopic 
insight into the driving forces for preferential adsorption or 
separation, for example by seeing how the adsorption of one 
component of the mixture enhances or suppresses the others, or 
by looking at guest-guest interactions. 

But probably the greatest interest in the modeling of coadsorp¬ 
tion through GCMC is in the computational prediction of fluid 
separation properties, which can be difficult and time-consuming 
to measure experimentally. This is particularly true of coadsorp¬ 
tion because of the increased dimensionality of the problem, 
with a larger number of control parameters. Even for a binary 
mixture, coadsorption properties are a function of temperature, 
pressure and mixture composition, thus requiring tedious ex¬ 
perimental work to map out in large ranges of each parameter. 
Even more so for ternary and more complex mixtures, with the 
number of experimental measurements for different composi¬ 
tions increasing exponentially with the number of components. 
There is thus great incentive in using molecular simulation of 
coadsorption to map out the separation properties of MOFs, find¬ 
ing optimal working conditions for applications, or providing 
thermodynamic information as input for process engineering and 
optimization. The process of a typically computational study 
of gas coadsorption through GCMC simulations of mixtures is 
schematized in Figure Force fields (an input of the computa¬ 
tion) are validated for each separate gas by comparison between 
experimental pure-component isotherms and GCMC simulations. 
If the agreement is not satisfactory, the MOF-adsorbate force 
fields need to be optimized: either adjusted empirically to the 
experimental data, or based on ab initio reference data. Once the 
force fields have been validated for pure components, GCMC 
simulations of the mixture at various temperature, pressure and 
composition are performed. Analyzing these, the results ob¬ 
tained include coadsorption isotherms, selectivities, heats of 
(co)adsorption, plot of adsorbed densities for each component, 
etc. A rather typical example of this type of study is shown in 
Figure for the case of coadsorption of xenon and krypton in 
MOF CPO-27-Ni. 12361 


The other alternative to the modeling of coadsorption in 
MOFs, especially for gas mixtures, is the use of analytical 
mixture models such as the Ideal Adsorbed Solution The¬ 
ory HAST). 12371 or more complex variants such as Real Ad¬ 
sorbed Solution Theory (RAST1 B238[ 12391 or Vacancy So¬ 
lution Theory II24Q[ I241II . While those have been shown 
to present reasonable predictions among mixtures of small 
molecules. 12421 [2431 1244II they can deviate from GCMC pre¬ 
dictions and experiments in the cases of highly-competitive ad¬ 
sorption, tight packing of guest molecules (shape selectivity) or 
MOF structures with marked chemical heterogeneity. ifTslI^ 
Thus analytical mixture theories, such as I AST, should only 
be used as a rough estimate of separation properties, before 
more time-consuming molecular simulations or experiments are 
performed. 

4.4. Beyond classical force fields: open metal sites and 
chemisorption 

The molecular simulation of adsorption by Grand Canonical 
Monte Carlo relying on classical force fields is a very pow¬ 
erful tool, but it is limited by the accuracy of a given force 
field to describe host-guest interactions in the system. There 
are several cases where the classical approximations are not a 
realistic description of the intermolecular interactions, in par¬ 
ticular in the case of chemisorption (i.e., formation of bond(s) 
between host and guest). Within the field of adsorption in metal- 
organic frameworks, the most common feature which classical 
force fields typically fail to describe is the short-range interac¬ 
tion of guest molecules with metal centers, i.e. the interaction 
with undercoordinated metal sites (or coordinatively unsaturated 
metal sites, cus). This generic limitation has been evidenced 
on several pairs of adsorbents and guest molecules, includ¬ 
ing hydrogen. 12461 methane. 12471 carbon dioxide. l248[ 12491 
propane and propylene. 12501 etc. 

In such cases, in order to accurately describe specific host- 
guest interactions (including guest interactions with coordina¬ 
tively unsaturated metal sites), one needs to use ah initio quan¬ 
tum chemistry methods, either in the Density Functional The¬ 
ory (DFT) approach, using post-Hartree-Fock methods, or multi¬ 
reference methods. The latter two families of methods can be 
very accurate and require no ad hoc parameters: they are thus 
very powerful to study novel systems for which classical force 
fields are not available, or where the adsorption mechanism is 
unknown, as well as very specific interactions. On this extensive 
literature, we refer the reader to the recent review of Odoh et 
al.d on quantum-chemical characterization of MOF proper¬ 
ties, which includes an extensive section on the investigation of 
adsorption properties. 

A rather representative example of this is the calculation of 
CO 2 adsorption in metal-organic frameworks of the CPO-27 
family with different metals, using a combination of DFT and 
post-Hartree-Fock methods. ESD Yu et al. have reported very 
accurate calculations of the CO 2 binding energies on open metal 
sites in CPO-27-M, with M = Mg, Mn, Fe, Co, Ni, Cu, or Zn. 
They first performed geometry optimizations of the periodic 
structures (atomic positions and lattice parameters) at the DFT 
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Figure 18: Flowchart of a typical computational study of gas mixture coadsorption in metal-organic frameworks through GCMC simulations. 
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Figure 19: Results from a computational co-adsorption study of xenon/krypton mixtures in metal-organic framework CPO-27-Ni, by Gurdal et al. (236) From left to 
right: single component and equimolar Xe/Kr mixture adsorption isotherms; adsorption selectivities at 1 and 10 bar, as a function of mixture composition; Xe/Kr 
adsorption and permeation selectivities as a function of pressure; performance of CPO-27-Ni for Xe/Kr separation compared to other MOFs. Adapted with permission 
from ref. ma. Copyright 2012 American Chemical Society. 


level with generalized gradient approximation (GGA) exchange- 
correlation functional. Adsorption geometries were determined 
by geometry optimizations on representative clusters of the ma¬ 
terials (see Figure]^, again at the DFT level. Accurate binding 
energies were then calculated using these geometries by single 
point energies at the MP2 level of theory. They were then fur¬ 
ther corrected with a QM/MM approach, where the QM energies 
were from the MP2-level calculations and the MM energies were 
calculated with the Dreiding and UFF force helds. 

In contrast with this high-accuracy methodology, studies of 
adsorption in MOFs with coordinatively unsaturated sites can 
also be performed entirely at the DFT level. These are com¬ 
putationally much cheaper and have been widely used in the 
literature. They need, however, to be carefully benchmarked 
against experimental data or highly accurate reference calcula¬ 
tions, as the choice of DFT exchange-correlation function and 
dispersion correction scheme can heavily influence the prop¬ 
erties calculated, such as adsorption energies. There appears 
to be no “one size fits all” choice of methodology in this area, 
and optimal exchange-correlation functionals for each adsor- 
bate/MOF pair need to be found by comparing the performance 
of various functionals on smaller cluster models using post- 
Hartree-Fock or multi-reference calculations. 112531 [25411252II 
This is shown for example on Figure a comprehensive plot 


of binding energies for several small molecules (methane, car¬ 
bon dioxide, propane, and propene) in HKUST-1, evaluated with 
various DFT exchange-correlation functionals and compared to 
DFT/CC results.ll252l 

Energy minimization calculations, as described above, only 
describe well-defined adsorption sites of a single molecule, as 
well as the corresponding geometries and low-coverage adsorp¬ 
tion enthalpies (i.e., adsorption enthalpies in the limit of zero 
gas pressure). They do not account for entropy and guest-guest 
interactions, and thus cannot treat pore filling or packing effects 
of the molecules inside the pores. As a consequence, they cannot 
describe the adsorption at high uptake or where entropy plays 
a big role. Therefore, there has been a very large effort in the 
literature to use the insight and data gained from ab initio calcu¬ 
lations (adsorption sites, binding geometries, energies, forces, 
etc.) in order to parameterize host-guest force fields. This 
can take the form of reoptimizing existing force fields, adjust¬ 
ing some of the terms to match the quantum chemistry data: 
typically, strengthen the metal-guest potentials to more accu¬ 
rately describe the binding of molecules to the open metal sites 
of the framework. It can mean choosing an entirely altogether 
functional form for the metal-adsorbate interactions, or even a 
numerical description based on the energy profile as a function 
of adsorbate-metal center distance. 
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Figure 20: Calculations of the CO 2 binding energies on open metal sites in CPO-27-M, with M = Mg, Mn, Fe, Co, Ni, Cu, or Zn. Adapted from ref. (SSD with 
permission from The Royal Society of Chemistry. 



Figure 22: Ab initio parameterization of a force field describing the adsorption of CO 2 in CPO-27-Mg. Adapted with permission from ref. 12491 . Copyright 2012 
American Chemical Society. 


There are several examples of force held optimization for 
open metal site MOFs in the literature, including the cases of: 
hydrogen in HKUST- 1: 112461 water in copper-based MOFs: ll255l 
carbon dioxide and water in Mg-MQF-74 I256I and Zn-MOF- 
74: 1125711 CO 2 adsorption in Fe^tdobdc): 1125811 etc. We detail 
here by way of illustration the procedure followed by Chen 
et al. ll249l for the ab initio parameterization of a force held 
describing the adsorption of CO 2 in CPO-27-Mg. Starting from 
a Buckhingham-type potential (the Carra-Konowalow potential), 
the authors calibrated a MMSV (Morse-Morse-spline-van der 
Waals) piecewise interaction potential for interactions with the 
MOF’s coordinatively unsaturated metal sites. The optimization 
of the MMSV potential was performed against reference ab 
initio data obtained with a double-hybrid density functional with 
empirical dispersion correction, B2PLYP-D2. II259[ [26011 The 
optimization itself was performed using a multiobjective genetic 
algorithm, (26D and the good agreement between energy curves 
and ab initio data is evident (see Figure [22|). Chen et al. then 
used the optimized force held to perform GCMC simulations of 
adsorption isotherms, and showed vastly improved agreement 
compared to the standard force helds such as UFF. 

Another solution that has been proposed is to precompute the 
MOF-guest interactions on a hne mesh of points within the pores 
of the MOF at the quantum chemical level. Using single-point 
quantum chemistry calculations, typically at the DFT level, the 
potential energy surface for a single guest molecule can be ob¬ 
tained, and the combined with guest-guest classical interatomic 


potentials in a series of GCMC simulations. This approach 
has been demonstrated by Chen et al. on the study of methane 
adsorption on HKUST-1 (also known as Cu 3 (btc) 2 ), a MOF 
with coordinatively unsaturated metal sites. 12471 This method 
combines a DFT level of description of the host-guest inter¬ 
actions with a full sampling of the phase space of the system, 
thus accounting for entropy. It is, however, only possible for 
atomic guests (including rare gases) or small spherical molecules 
(such as methane): for nonspherical molecules, the additional ro¬ 
tational degrees of freedom make the approach computationally 
prohibitive. 

4.5. Adsorption in flexible MOFs 

Although Monte Carlo simulations in the Grand Canonical 
ensemble, as presented at the beginning of this section, are 
considered the gold standard in the simulation of adsorption in 
nanoporous materials, they rely on a very strong assumption: 
that the host material is rigid. This approximation is very reason¬ 
able when it comes to the thermodynamics of adsorption in most 
porous inorganic materials, such as zeolites, where framework 
flexibility is limited. However, because metal-organic frame¬ 
works are based on weaker bonds and interactions (coordinative 
bonds, n-n stacking, hydrogen bonds, etc.), that are responsible 
for their intrinsic structural flexibility. The organic-inorganic 
connections therefore allow underconstrained structural linkages 
that are responsible for soft mechanical properties, and organic 
linkers with side chains allow for local dynamics of the host 
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Figure 21: Errors in the interaction energies calculated for methane, carbon 
dioxide, propane, and propene at CUS (open metal site), WIN (cage window), 
and CTR (cage center) sites of HKUST-1 with respect to the DFT/CC reference 
level of theory. Adapted with permission from ref. 12521 . Copyright 2015 
American Chemical Society. 


framework. This flexibility of MOFs can be triggered upon ad¬ 
sorption, leading to large-scale structural changes in the MOF 
framework and subsequent alteration of its physical and chemi¬ 
cal properties, including adsorption itself. Such flexible MOFs, 
sometimes called Soft Porous Crystals, OlTl are in growing 
number. 1262112631 ri48l There is therefore a need for molecular 
simulations methods describing the interplay between adsorption 
and host flexibility, going beyond the GCMC method. When the 
flexibility occurs in the form of a clear transition between well- 
defined and identified crystallographic structures, a simple “two 
state” solution can be used, studying the adsorption in both rigid 
structures through GCMC. 112641126511 However, more complex 
systems require a more direct approach, accounting directly for 
their flexibility in molecular simulations. 

From a thermodynamics point of view, the adsorption of 
molecular fluids inside deformable hosts is most appropriately 
described in the osmotic thermodynamic ensemble, where the 
control parameters are the number of molecules of the host 


framework A^host^ the chemical potential of the adsorbed fluid 
/iads. the mechanical constraint G exerted on the system and the 
temperature T. (266) Direct molecular Monte Carlo simulations 
of adsorption in this open ensemble are possible, where the 
Monte Carlo moves typically used in GCMC are supplemented 
with “volume change” moves (see Figure [T7]). However, the 
efficiency of these simulations are greatly improved by perform¬ 
ing simulations in a hybrid MC/MD setup (also called Hybrid 
Monte Carlo, or HMC).||267||268ll263 In HMC, short molecu¬ 
lar dynamics trajectories are considered as Monte Carlo moves, 
allowing to better sample the host framework’s flexibility by 
following its collective motions. 

Maybe one of the most striking example of osmotic Monte 
Carlo simulations is that reported by Ghoufi et al. on the CO 2 
adsorption-induced breathing of MIL-53(Cr). ll270[ IZTlI This 
study showed that HMC in the osmotic ensemble was able to 
describe both thermal, mechanical and adsorption-induced struc¬ 
tural transitions between the two phases of MIL-53(Cr), based 
on an ad hoc force field describing the material (as depicted 
in Figure [23] ). This study also demonstrated the main issue 
with direct HMC simulation, namely the widely hysteretic na¬ 
ture of the transitions and the difficulties in overcoming free 
energy barriers associated with the transition and determining 
the thermodynamic equilibrium of the system. 

Another way to perform molecular simulation in the osmotic 
ensemble, avoiding the convergence issues of the direct Hy¬ 
brid Monte Carlo method, is to rely on free energy methods to 
calculate the osmotic potential of the system for each possible 
value of chemical potential. These methods require the use of 
non-Boltzmann sampling in (guest loading, volume) parameter 
space in order to fully characterize the adsorption thermody¬ 
namics as well as the material’s response to adsorption. In 
addition, they provide information on both the thermodynamic 
equilibrium as well as all metastable states of the system. Three 
different variants have been proposed: via thermodynamic inte¬ 
gration based on GCMC simulations as performed by Watanabe 
et al.. ll272[|Z7^ through the Wang-Landau algorithm as done 
by Bousquet et al.. 11274 [|Z75]| or similarly with the Transition- 
Matrix Monte Carlo sampling as proposed by Shen et al. 12761 
All three methods were demonstrated on model systems, and 
although they appear promising they have not yet been applied 
to atomistically-detailed molecular frameworks. 

Finally, a third way to model adsorption in Soft Porous Crys¬ 
tals is the use of thermodynamics-based analytical models (see 
Ref. m for a review on this topic). This is achieved by writ¬ 
ing down the equations for adsorption in the osmotic ensem¬ 
ble and introducing simple approximate expressions for cer¬ 
tain quantities, such as describing the adsorption isotherms 
can be modeled by classical equations such as Langmuir or 
Langmuir-Freundlich, which can be either fitted from experi¬ 
mental data ll266[ I277II computed from GCMC simulations of 
rigid frameworks. l278[[Z79l 12801 or calculated through param¬ 
eterized equations of state. 1128111 This approach was first used 
to study stepped adsorption isotherms in bistable materials of 
the gate opening and breathing families. II266I I282ll and later 
extended to take into account the influence of temperature 1128 3 II 
and mechanical pressure. GSD 
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Figure 23: Left: the two metastable phases of “breathing MOF” MIL-53(Cr): large-pore (top) and narrow-pore (down). Right: structural transitions observed by 
Hybrid Monte Carlo molecular simulations, under stimulation by temperature, mechanical pressure, and carbon dioxide adsorption. Adapted with permission from 
refs. 127H and 12701 . Copyright 2010, 2012 American Chemical Society. 



Figure 24: From a set of experimental adsorption isotherms in a flexible MOF, here for xenon in the “breathing” framework MIL-53(A1), analytical thermodynamic 
models allow the construction of a full (temperature, gas pressure) phase diagram. Reproduced from ref. (283). Copyright 2009 Wiley-VCH Verlag GmbH & Co. 
KGaA, Weinheim. 


Such analytical models not only help rationalized the behav¬ 
iors observed experimentally, and shed light into their key ther¬ 
modynamic factors and driving forces, but they can also have 
predictive value. This is, in particular, the case of the OFAST 
model l284[ 12851 (osmotic framework adsorption solution the¬ 
ory), extending the widely-used I AST coadsorption model to the 
osmotic ensemble. Based solely on experimental pure compo¬ 
nent adsorption isotherms, the OFAST model is able to predict 
coadsorption in flexible MOFs, and has been fully validated by 
comparison against experimental data, e.g. on the coadsorption 
of CO 2 /CH 4 mixtures in MIL-53(AD. 12861 

4.6. Comparing simulations and experiments 

This short review on the computational modeling of adsorp¬ 
tion in MOFs cannot be complete without a few words on the 
comparison between computational and experimental data — 


and the pitfalls thereof. It is always desirable to validate a simu¬ 
lation methodology by comparing its results to available exper¬ 
imental data, when possible. However, the ample literature on 
MOF adsorption shows that simulation results and experimental 
data very often differ quantitatively 12 87 1 (as was already noted 
in the early work in the area ll288ll ). Moreover, the differences 
can be in some cases rather large, without the simulation being 
necessarily at fault: there is also large variability in experimen¬ 
tal measurements of adsorption data, and primarily adsorption 
isotherms. The possible causes for this variability are many: 

1. Dependence on the measurement technique used: volu¬ 
metric vs. gravimetric, equilibration times, etc. 12891 [2901 
Experimental isotherms may in some cases not truly rep¬ 
resent the thermodynamic equilibrium, especially when 
adsorption kinetics is slow (e.g., bulky molecules in small 
pore MOFs, or occurrence of structural transitions). 
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2. Impact of the activation procedure and dependence on the 
history of the sample. Residual amounts of solvent, organic 
linker or templates in the pores of a MOF can severely 
diminish its accessible pore volume and specihc surface 
area. 

3. Influence of the textural properties, with key parameters 
such as crystal size distribution, nature of the external sur¬ 
faces, crystal shape, etc. Adsorption on MOF samples of 
nanoscopic sizes occurs both inside the pores and on the 
external surface, (2791 or at crystallographic line or plane 
defects. Molecular simulations almost always assume a 
bulk behavior (i.e. an infinite crystal), while crystal sizes 
in MOFs vary widely depending on material, synthesis 
conditions, etc. 

4. Presence of defects, in particular missing organic linkers 
in the crystalline structure. Depending on synthesis con¬ 
ditions, some MOFs can present large amounts of defects 
within their crystalline structure. These defects can have a 
dramatic impact on adsorption properties, typically increas¬ 
ing the pore volume and pore sizes. (831 1 1 261129 1 1 12921 but 
also affecting the host-guest interactions and the chemical 
nature of the internal surface of the material. 12931 More¬ 
over, if the defects are organized rather than randomly 
dispersed throughout the MOF structure (e.g., correlated 
disorder in UiO-66 12941 ). the effects on adsorption can be 
expected to be even more important. 

For all these reasons, it is highly recommended to perform 
a full characterization of the pore volume of MOFs, compar¬ 
ing their geometrical properties to the experimentally measured 
surface areas, pore volume, pore size distribution, etc. 12951 l69l 
This is a necessary prerequisite to confirm that the simulations, 
performed on a perfect crystal structure, can match the experi¬ 
mental system and that any comparison of adsorption isotherms 
is meaningful. In some cases where the geometrical and exper¬ 
imental surface areas differ, it is possible to adjust the simula¬ 
tion results by scaling the adsorbed quantities by an empirical 
factor. 1216[ 12961 The reasoning behind this is to account for 
changes in pore volumes through blocked pores or missing link¬ 
ers, which are probably the most common issues. The scaling 
factor is thus dependent on the material studied, the synthesis 
and activation conditions, the exact form and textural proper¬ 
ties of the sample used experimentally... Scaling factors used 
in the literature typically vary between 0.7 and 1. The use of 
such a scaling procedure allows one to compare or extrapolate 
computational data for different guests in the same adsorbent 
material. However, this “quick fix” cannot truly replace a better 
understanding of the origins of the discrepancy. 

In conclusion, caution should be taken in comparing simulated 
adsorption results with experimental ones, and in particular 
when adjusting computational parameters (such as force field 
parameters) to fit certain experiments: one should do so not 
on a few isolated adsorption isotherms, but only based on a 
comprehensive experimental dataset (adsorption isotherms in 
wide temperature and pressure ranges, heats of adsorption, etc.). 


5. Perspectives 

Ending this introductory review of computational studies of 
metal-organic frameworks and their properties, we highlight 
some of the very recent work and remaining open questions that 
seem, from our perspective, to be important challenges for the 
development of the MOF field. 

5.1. Modeling of defects and disorder 

The occurrence of defects and their correlations have long 
been recognized to play a central role in the physical and chemi¬ 
cal properties of materials. The study of defects in crystalline 
compounds forms an important part of solid-state chemistry and 
physics, as can be seen by the large body of literature concern¬ 
ing dense and porous inorganic materials. 12971 12981 Yet, the 
importance of the role of defects and disorder in metal-organic 
frameworks is only starting to emerge, and is still rarely stud¬ 
ied in the existing literature, both experimental and theoretical. 
Yet, there is evidence that defects and disorder play an impor¬ 
tant role in the function of several existing MOFs. The poster 
child for this is the UiO-66(Zr or Hf) family of structures, of 
high interest because of their thermal, mechanical and chemical 
stability. Studies in the past three years have shown that UiO- 
66 materials can contain significant amount of missing-linker 
defects. 1 126112991 which have a crucial impact on the adsorption 
and catalytic properties of the material. II 3 0011 The concentration 
of these defects can be tuned during synthesis; C261 they do not 
occur randomly in the structure but are correlated and form nano- 
scaled domains of well-defined {reo) topology. 12941 Moreover, 
the introduction of controlled heterogeneity in MOFs, without 
loss of its ordered structure, can lead to the creation of more com¬ 
plex structures and pore environments. BO 11I302II introducing ad¬ 
ditional functions into known topologies and structures. (41 13031 
A recent review on the topic can be found in Ref. B04II . 

In stark contrast to the few examples given above, most MOF 
structures reported in the literature feature a crystallographic 
structure and basic characterization of porosity, as well as some 
macroscopic study of their function (adsorption, catalytic activ¬ 
ity, etc.). Because computational chemistry techniques are based 
on periodic representations of the crystalline structures, they 
model perfect materials with no defects and no disorder. A col¬ 
league summarized this first-order approach, tongue-in-cheek, 
as follows: for the most part, we are dealing with defects by 
‘'ignoring them, but invoking them as reason for any difference 
between modeling and experiments’’ 

Yet, some recent studies have used quantum chemistry and 
molecular simulation techniques to understand and predict the 
impact of defects on the properties of MOFs. At the quan¬ 
tum chemical level, a good example is the study by Chizallet 
et al. using DFT calculations of the catalysis of transesterifi¬ 
cation in ZIP-8, on acido-basic sites located at defects or at 
the external surface of the material. B06I Another one is the 
influence of defect concentration on the structural, mechani¬ 
cal and thermal properties of UiQ-66(Hf). 12941 giving micro¬ 
scopic insight into the occurrence of defect-dependent thermal 
densification and colossal negative thermal expansion (NTE) 
measured experimentally. B07II At the classical level, Ghosh et 
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al. have studied the impact of defects on adsorption properties 
of water in UiO-66(Zr) through Grand Canonical Monte Carlo 
simulations. 12931 

Though there have been a few studies (listed above) on the 
impact of defects on MOF properties, it is worth noting that the 
bigger question has not been addressed so far from the theoretical 
point of view: what causes certain MOFs to present defects, 
and when and why are these defects correlated? There is 
experimental evidence that the occurrence of defects in MOFs 
can be random in some cases, and correlated in others. For 
example, structures obtained by linker-exchange BOSl 13091 show 
no correlation between the substituted linker positions, i.e. the 
resulting partially solvent-exchange MOF is a solid solution of 
the two parent MOFs, and features a random arrangement of 
linkers. 13101 On the other hand, missing-linker defects in UiO- 
66(Hf) exhibit correlated disorder, with the occurrence of nano- 
scaled domains of weh-dehned topology. 12941 Modeling studies 
have, so far, been unable to address this kind of very fundamental 
questions. Nor have there been any computational studies on 
noncrystalline disordered MOF structures, i.e. molten, glassy, 
and amorphous MQFs. l[T^lmll3T2ll3T3ll 

5.2. Databases for high-throughput screening 

One of the areas that has seen fast-paced development in the 
last few years is that of high-throughput computational screening 
of porous materials in general, and metal-organic frameworks 
in particular. As described in some more detail in Section [23} 
this is the result of the conjunction of a few factors: (i) a large 
number of know MOF structures, (ii) the availability of both 
computational power and advanced modeling techniques, and 
(iii) a political incentive to speed up the discovery of novel 
materials. This has lead to the development of two different types 
of large-scale databases of MOF structures: hypothetical MOF 
structures generated by combinatorial approaches (of the order 
of 100k structures); (501 and “computation-ready” structures 
derived from experimental crystallographic data (of the order 
of 5k structures). ED The availability of these databases offers 
great opportunities for high-throughput computational screening 
of materials for specific applications, as has already been done 
for adsorption of hydrogen, ESED methane. (501 [65l carbon 
dioxide, E3 noble gases, (631 [64l etc. But this also raises some 
novel questions and challenges about how best to exploit these 
databases for materials discovery. 

Focusing first on databases of hypothetical structures, it might 
be here interesting to draw some parallels with the existing 
databases of hypothetical zeolitic structures. (5^ l55l which have 
existed for nearly ten years now and are of similar size (up to two 
million structures). There also, like for MOFs, the overwhelm¬ 
ing majority of screening studies have focused on adsorption 
properties, relying both on computationally-cheap geometrical 
characterization (pore size distributions, pore diameters, surface 
area, pore volume; see Section [23] ) and GCMC calculations 
on candidate structures selected on geometric criteria. This 
approach has yielded valuable insight into the relationship be¬ 
tween geometric properties and adsorption performances. (50l as 
well as the intrinsic limits of these materials. (651 H is, however, 
rather uncertain how this approach can be extended to the study 


of other properties, either for MOFs or zeolites: simple and 
“cheap” descriptors are hard to identify for other functions 

such as catalytic activity, thermal and mechanical properties, 
luminescence, fiexibility, ... 

Furthermore, given the large number of structures available, 
the use of a few descriptors and their correlation to key prop¬ 
erties is typically an overdetermined problem. Even simply 
plotting the data generated can become difficult for multidimen¬ 
sional data sets: two-dimensional correlation plots, for example, 
are limited to testing hypotheses formulated a priori. There is 
thus a real need to use more sophisticated data mining tools 
and to move away from predetermined descriptors. Some 
steps have been taken along this path in the use of quantita¬ 
tive structure-property relationship (QSPR) backed by machine 
learning algorithms, for example by Fernandez et al. to iden¬ 
tify high-performing MOFs for carbon dioxide capture. Eia 
Thornton et al. used a combination of molecular simulation and 
machine-learning techniques to identify candidate zeolites for 
catalytic reduction of carbon dioxide, from a sample of 300 thou¬ 
sand zeolite structures (hypothetical and experimental). (am At 
a much smaller scale, Yildiz et al. have demonstrated the poten¬ 
tial of Artificial Neural Networks (ANN) to predict adsorption 
properties on the case of hydrogen gas storage in thirteen MOFs 
with high surface areas. (3T61 

Another issue that needs to be addressed is that of the ex¬ 
perimental feasibility of the hypothetical MOFs. Hypothetical 
zeolite databases have existed for quite some time, yet from their 
systematic exploration no real structure has been synthesized and 
tested for practical applications, and the question of why so few 
zeolites are observed while there are so many hypothetical frame¬ 
works remains open. l3171l3181l319l Though MOFs seem to far 
a little better here, with several examples of materials synthe¬ 
sized after they have been selected by computational design, Ei 
the important questions of experimental feasibility and compu¬ 
tational rational design remain. Among the currently-listed 
hypothetical MOF structures, which are experimentally acces¬ 
sible to synthesis and present sufficient mechanical, thermal 
and chemical stability for practical applications? Can we 
use the tools of theoretical chemistry to help guide the synthesis 
of novel MOFs identified for their properties? Can we help pre¬ 
dict the conditions (reagents, solvents, temperature, etc.) under 
which a given MOF may be obtained experimentally? 

Finally, databases built from experimental structures, such 
as CoRE M0E. (6TI| raise important questions as to the cura- 
tion of the database and the relevance of the “cleaned-up” 
structures to the properties of the original experimental ma¬ 
terials. The generation of computation-ready structures from 
crystallographic data involves an automated procedure to re¬ 
move solvent and disordered guests. This procedure has severe 
limits. Eirst, it assumes that all coordinated solvent molecules 
can be removed from the structure: this is true in some cases 
(water molecules bound to HKUST-1 metal centers, for exam¬ 
ple) but not in general. Indeed, there exist some MOEs with 
included solvent where the very stability of the framework is 
crucially dependent on the presence of the solvent molecules 
(e.g. water, methanol, etc.), and which cannot be activated with¬ 
out loss of crystallinity. II320[ [3211 1 16611 Secondly, it assumes 
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that there is no relaxation upon solvent removal, which is a 
rather crude approximation, especially in the case of soft porous 
crystals. Nonetheless, this large-scale database is the hrst of its 
kind and is likely to prove useful for understanding fundamental 
principles of MOFs and a basis for computational selection of 
structures for targeted applications. 


5.3. Stimuli-responsive MOFs 

With the large number of new MOF structures synthesized and 
characterized every year, one of the empirical patterns appearing 
is the common occurrence of flexibility of these framework 
materials. There is a rapidly increasing number of framework 
structures whose flexibility manifests in the form of large-scale 
structural transformations induced by external stimulation of 
physical or chemical nature: changes in temperature, mechanical 
constraints, guest adsorption, light exposure, etc. A number of 
different terms have been used to refer to this behavior, including 
smart materials, soft porous crystals, ClTj dynamic frameworks, 
flexible frameworks, and stimuli-responsive materials. (Ml The 
realization, over the past ten years, of the prevalence of these 
flexible materials and their potential for applications has lead to 
the emergence of an entire subfield of computational methods 
dedicated to their modeling. The existing literature on this 
topic has been described in Sections |4.5| (adsorption-induced 


structural transitions) and Section [3T] (mechanical properties). 
This is, however, an area in fast development and with many 
challenges still left open. 

The first is the need for the systematic development of high- 
accuracy force fields for highly fiexible materials. Force 
fields for MOFs with high degree of flexibility, such as MIL-53 
or ZIFs, require a good accuracy in the description of the in¬ 
tramolecular interactions of the framework, and in particular the 
low-frequency phonons modes. For this, two categories of force 
fields have been used so far. The first, and dominant category, 
is the use of transferable intramolecular force fields for organic 
molecules (e.g., from the AMBER parameter database 13221 ) 
along with hand-tuned metal-organic bending and torsion po¬ 
tentials. The later are then optimized to reproduce known struc¬ 
tural properties and vibration frequencies, 11323 [ I324[ 13251 or 
experimental data such as adsorption isotherms and structural 
transitions. I326[|327l 13281 This approach is not entirely gener- 
alizable, and the quality of the resulting force field is only vali¬ 
dated on limited experimental data and cannot be systematically 
improved, as the optimization problem is typically undercon¬ 
strained (too many parameters fitted on relatively little data). On 
the other hand, force fields derived from ab initio data (as de¬ 
scribed in Section [44| ) show clear, but that methodology is rather 
hard to apply to flexible materials, especially when it comes to 
the intramolecular terms. Still, some recent progress has been 
made in this area, giving hope that generic methodologies for 
ab initio force fields of flexible MOFs can be a reality in the 
future. 123II 

The second area of development we want to highlight here 
is the recent trend towards better insight into the microscopic 
mechanisms of stimuli-responsiveness, with particular focus on 
space-resolved, time-resolved, and in operando experimental 
measurements. Those are crucial to provide a better fundamental 


understanding of the fundamental nature of the transformations 
triggered by complex stimuli and have practical consequences 
for the design of novel materials in working conditions. Yet, 
relatively little theoretical and computational effort has been 
spent on this, most probably due to the very difficult nature 
of the issue. The questions that need to be answered include: 
how do stimuli-induced transformations occur and propa¬ 
gate at the scale of the crystal? What is their kinetics and 
dependence on the history of the material? What determines 
the possible metastable states of the system? How do crystal 
size, shape, and textural properties affect their physical and 
chemical properties, and ultimately their responsivity? 

Some of the recent studies, both experimental and theoret¬ 
ical, have started to address this issue. On the experimental 
side, one striking such example was the recent observation of 
an transient state in the adsorption-induced structural transition 
of pillared MOF Zn 9 (ndc) 9 (dabco), 13291 by using synchrotron 
grazing incidence diffraction measurements to determine sepa¬ 
rately the structures of a crystal’s bulk and surface. II330I Kondo 
et al. demonstrated that, upon adsorption of bulky slow-diffusing 
molecules, the MOF featured a heterostructure with a guest- 
induced sheared phase at the surface coexisting with an un¬ 
perturbed MOF structure in the core of the crystal. On the 
theoretical point of view, we presented some time ago a multi¬ 
scale physical mechanism and a stochastic model of breathing 
transitions. 1133II which is to our knowledge still the only ex¬ 
ample of such “upscaling”, trying to address computationally 
the dynamics of structural transitions at the level of the crystal 
itself. Based on a simple Hamiltonian that describes the physics 
of host-host and host-guest interactions, we showed how the 
behavior of unit cells is linked to the transition mechanism at the 
crystal level through three key physical parameters: the transi¬ 
tion energy barrier, the cell-cell elastic coupling, and the system 
size. (3321 

Along a similar line, but on the different topic of MOF crys¬ 
tal growth, Yoneya et al. have recently used coarse-grained 
molecular dynamics to investigate the MOF self-assembly and 
the influence of metal-ligand coordination in this process. 13331 
This is a first step toward a better understanding of MOF synthe¬ 
sis and growth at the microscopic scale, a topic that has proven 
crucial but difficult to address so far, in MOFs as well as in other 
porous materials. 1334113351 
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HMC: 
HPC: 
lAST: 
IRMOF: 
MC: 
MD: 
MMSV: 
MOF: 
NFC: 
NTE: 
OFAST: 

PSD: 

QSPR: 

RAST: 

RCSR: 

RESP: 

SBU: 

SPG: 

TDDFT: 

UA: 
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Cambridge Structural Database 
Covalent Organic Framework 
Computation-Ready Experimental MOF 
Density Functional Theory 
Grand Canonical Monte Carlo 
Hybrid Monte Carlo 
High-Performance Computing 
Ideal Adsorbed Solution Theory 
Iso-Reticular Metal-Organic Framework 
Monte Carlo 
Molecular Dynamics 
Morse-Morse-spline-van der Waals 
Metal-Organic Framework 
Negative Finear Compressibility 
Negative Thermal Expansion 
Osmotic Framework Adsorption 
Solution Theory 
Pore Size Distribution 

Quantitative Structure-Property Relationship 
Real Adsorbed Solution Theory 
Reticular Chemistry Structure Resource 
Restrained Electrostatic Potential 
Secondary Building Unit 
Soft Porous Crystal 

Time-Dependent Density Functional Theory 

United Atom 

Universal Force Field 

Vacancy Solution Theory 

Zinc Alkyl Gate 

Zeolitic Imidazolate Framework 

Zeolite-like Metal-Organic Framework 
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